#!/usr/bin/env bash

###check for parameters' number
if [ $# -lt 4 ];then
echo -e "error : liftover.sh mergefile databasefile date_of_freqfiles file1 file2 ... \n file1 = liftover_wellconverted_dataname1.txt, file2 = liftover_wellconverted_dataname2.txt etc..." 1>&2
exit 1
fi

#test=`find db*`
#if [ $? -ne 0 ]; then
#echo "ce fichier n'est pas une image .jpg"
#else
#/usr/bin/convert $f -resize $geom ~/Images/thumbnails/$(basename $f .jpg)_$geom.png
#fi
#find db* > /dev/null
#if[ $? -ne 0 ];then

datefile=`date +%e%m%Y` # date ddmmyyyy

### mergefile
mergefile=$1
shift
cat $mergefile | sort -t' ' -k1,1 > merge_tmp1.txt

### databasefile

db=$1
cat $db |  awk -F" " '{print $1"-"$3,$1,$2,$4,$5}' | sort -t' ' -k1,1 > "db_"$datefile".txt" #sort 1000G database file by chromosome and position
shift

### create output directory if it doesn't exist yet
if [ ! -d "./output_liftoversh" ]; then
	mkdir output_liftoversh
fi

# frequencies files date
datefreqfile=$1
shift

for f in $@;do

	echo "FILE : $f " 
	### datafiles wellconverted "rs" SNPs
	name=`basename $f .txt` # remove ".txt" suffix
	cat $f | sed 's/^chr//' | awk -F"	" '{print $1"-"$2,$4}' | sort -t' ' -k1,1 > "sort_"$name"_"$datefile".txt" # sort by chromosome and position
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/summary_"$name"_"$datefile".txt" # header
	join -t' ' -1 1 -2 1 "sort_"$name"_"$datefile".txt" "db_"$datefile".txt" | tr ' ' '	' >> "output_liftoversh/summary_"$name"_"$datefile".txt"  # join files by identical chromosome-position

	# select lines with rsID in 1000genomes database 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB" > "output_liftoversh/rs_"$name"_"$datefile".txt"  # header 
	cat "summary_"$name"_"$datefile".txt" | egrep "	[0-9]+	rs" >> "output_liftoversh/rs_"$name"_"$datefile".txt" 

	echo "see output_liftoversh/rs_"$name"_"$datefile".txt"

	# datafiles wellconverted "chr-pos" SNPs and frequencies' files
	if [ ! -f "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out"];then
		echo "for frequencies_dataname.bim_dateoffreqfile.out liftover file must be : liftover_wellconverted_dataname.txt and liftover_failed_dataname.txt"
	exit 2
	fi
	cat "frequencies"${name#liftover_wellconverted}".bim_"$datefreqfile".out" | sort -t'	' -k2,2 | awk -F"	" '{print $2,$6,$7,$8,$4,$5}' | tr '	' ' ' > tmpfreq.txt
	# select lines with chrposID in 1000genomes database 
	cat "summary_"$name"_"$datefile".txt" | egrep -v "	[0-9]+	rs" | sort -t'	' -k2,2 | tr '	' ' ' > "tmpchrpos_"$name"_"$datefile".txt" 
	echo "chr-position	snpIDdata	chr	snpIDdb		a0DB	a1DB	freqDB_a0 freqdata_a0	deltafreq	a0data	a1data" > "output_liftoversh/chrpos_"$name"_"$datefile".txt" #header of chrpos file
	join -1 2  -2 1  "tmpchrpos_"$name"_"$datefile".txt" tmpfreq.txt | awk -F" " '{print $2,$1,$3,$4,$5,$6,$7,$8,$9,$10,$11}' | tr ' ' '	' >> "output_liftoversh/chrpos_"$name"_"$datefile".txt"  

	echo "see output_liftoversh/chrpos_"$name"_"$datefile".txt"

	rm "sort_"$name"_"$datefile".txt" 

	### datafiles failed
	endfailedfile=`basename ${f#liftover_wellconverted} .txt`
	failedfile="liftover_failed"$endfailedfile
	cat $failedfile".txt" | sed '/#Deleted in new/d' | cut -d"	" -f4 | uniq > "output_liftoversh/exclude_"$failedfile"_"$datefile".txt"
	echo "see output_liftoversh/exclude_"$failedfile"_"$datefile".txt"

	### check RsMerge
	# delete "rs" prefix 
	cat "output_liftoversh/rs_"$name"_"$datefile".txt" | sed ':z;s/rs//;tz' > temporary.txt
	# search rsID 1000G in RsMerge
	cut -d"	" -f4 temporary.txt | sort > listotest.txt
	join -1 1 -2 1 listotest.txt merge_tmp1.txt  > "jointmp2_"$name"_"$datefile".txt"
	cat "jointmp2_"$name"_"$datefile".txt" | sort -t' ' -k2,2 > "jointmp3_"$name"_"$datefile".txt"
	cut -d"	" -f2 temporary.txt | sort > listotest2.txt
	join -1 1 -2 2 listotest2.txt "jointmp3_"$name"_"$datefile".txt" | tr ' ' '	' > "output_liftoversh/keep_"$name"_"$datefile".txt"
	echo -e "see output_liftoversh/keep_"$name"_"$datefile".txt \n" #-e allows \n interpretation

done

rm tmpfreq.txt
rm temporary.txt
rm listotest.txt
rm listotest2.txt
rm jointmp*

# command example : ./liftover2.sh ../../db/RsMergeArch.bcp ../../db/NEW_1000/EUR.chrfulldata.db 7122011 liftover_wellconverted_*
